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Abstract 


As the second of Earth’s Trojan asteroids, 2020 XL; is worthy of rendezvous and even sample return missions 
in many aspects. In this paper, a rendezvous mission to Earth’s second Trojan asteroid 2020 XL; is proposed. 
However, due to its high inclination and large eccentricity, direct impulsive transfer requires large amounts of 
fuel consumption. To address this challenge, we explore the benefits of electric propulsion and multi-gravity 
assist techniques for interplanetary missions. These two techniques are integrated in this mission design. The 
design of a low-thrust gravity-assist (LTGA) trajectory in multi-body dynamics is thoroughly investigated, 
which is a complex process. A comprehensive framework including three steps is presented here for 
optimization of LTGA trajectories in multi-body dynamics. The rendezvous mission to 2020 XLs is designed 
with this three-step approach. The most effective transfer sequence among the outcomes involves Earth— 
Venus—Earth-Venus-2020 XL;. Numerical results indicate that the combination of electric propulsion and 
multi-gravity assists can greatly reduce the fuel consumption, with fuel consumption of 9.03%, making it a 
highly favorable choice for this rendezvous mission. 
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1. Introduction 


It was proved in 1772 that there are two groups of small 
bodies which can stably share the orbit of a planet if they 
remain near the triangular point 60° ahead (L, point) or behind 
(Ls point) it in its orbit (Connors et al. 2011). The first 
discovery of these Trojan asteroids dates back to 1906, when 
Max Wolf found 588 Achilles (Nicholson 1961), orbiting 
around the triangular point 60° ahead of Jupiter along its orbit. 
Since then, more than 10,000 Trojans have been identified for 
Venus (De la Fuente Marcos & de la Fuente Marcos 2014), 
Mars (de La Fuente Marcos & de La Fuente Marcos 2013), 
Jupiter (Dvorak & Schwarz 2005), Uranus (de la Fuente 
Marcos & de la Fuente Marcos 2015) and Neptune (Sheppard 
& Trujillo 2006), with the majority being Jupiter Trojans. 
However, no Earth Trojan (ET) was identified until 2011. 
Asteroid 2010 TK; was observed in the Wide-field Infrared 
Survey Explorer (WISE) mission, and was found to be the first 
ET asteroid (Connors et al. 2011). As the current limits allow a 
population of perhaps 10° ETs with diameters greater than 
100 m (H ~ 23) (Wiegert et al. 2000), there are more ETs to be 
discovered besides 2010 TK7. However, due to the unfavorable 
viewing geometry of an object orbiting Earth-Sun’s L4 or Ls 
points as seen from Earth (Malhotra 2019), all the dedicated 
surveys aimed at detecting new ETs found nothing (Wiegert 
et al. 1997; Whiteley & Tholen 1998; Markwardt et al. 2020; 
Lifset et al. 2021). 


Asteroid 2020 XL; was discovered by the Panoramic Survey 
Telescope and Rapid Response System (Pan-STARRS) survey 
on 2020 December 12,’ and was considered as a potential 
candidate of an ET asteroid (de la Fuente Marcos & de la 
Fuente Marcos 2021). Subsequent investigations, conducted 
through follow-up studies, solidified 2020 XL5’s identity as the 
second ET asteroid (Hui et al. 2021; Santana-Ros et al. 2022). 
Both 2010 TK; and 2020 XL; are classified as transient ETs, 
with their orbital stability around L4 shown to be on the scale of 
thousands of years, far from the stability timescale of a 
theoretical primordial ET population (Santana-Ros et al. 2022). 
The inferred diameter for 2020 XL; is 1.18t898 km with an 
assumed albedo of 0.06 + 0.03 (Santana-Ros et al. 2022). This 
diameter is larger than that of the first ET asteroid 2010 TK , 
which was estimated to have a diameter of ~0.3 km (Connors 
et al. 2011). Because of some non-gravitational effects like the 
Yarkovsky effect, the orbital motion is related to its physical 
and geological properties. The determination of such properties 
becomes crucial to study the origin and orbital evolution of the 
asteroid. Moreover, because ETs share the orbit of Earth, their 
chemical composition can greatly reflect the space environment 
of Earth. Therefore, ETs are significant candidates for 
rendezvous and even sample return missions. 

However, due to the relatively high inclination and 
considerable eccentricity of both 2020 XL; and 2010 TK, 
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reaching the targets directly from low Earth orbit (LEO) using 
chemical propulsion (CP) proves to be a costly mission. The 
minimum total delta-v for the rendezvous mission with ballistic 
trajectory varies between 7.9 and 10.3 kms! for 2020 XLs, 
depending on the launch conditions. Meanwhile, for 2010 TK, 
this range is between 6 and 8.5 kms ' (Santana-Ros et al. 
2022). Thus, a rendezvous mission to 2020 XL; is much 
costlier than a mission to 2010 TK 7. Considering the delta-v 
budget, it is imperative to further optimize the trajectory for the 
rendezvous mission. For a mission to 2010 TK 7, there have 
been relevant studies about the optimization of transfer 
trajectory (Lei et al. 2017; Gao et al. 2019). The gravity-assist 
technique can be adopted to reduce the launch energy and 
rendezvous delta-v since it can effectively change the 
inclination of a spacecraft’s trajectory. In two-body dynamics, 
numerical results show that the fuel consumption of an 
impulsive trajectory with a Venus-Earth-Venus swingby 
sequence is 41% with flight time of 1694 days (Lei et al. 
2017). Moreover, electric propulsion (EP) is also an efficient 
technique to save the fuel consumption (Chen et al. 2018). EP 
was applied in the last leg of an impulsive transfer trajectory, 
i.e., from Venus to 2010 TK,, transforming the impulsive leg 
into a low-thrust leg (Lei et al. 2017). The application of EP in 
the last leg reduced the fuel consumption to 4.59% with flight 
time of 1684 days, where the initial mass was 800 kg and the 
maximal magnitude of thrust was 120 mN. This work aims for 
the design of a rendezvous mission to the new ET asteroid 
2020 XLs. In this work, since the orbit of 2020 XL; has a high 
inclination together with a large eccentricity, EP and multi- 
gravity assist techniques are combined for the optimization of 
transfer trajectory. 

There have been numerous studies about the optimization of 
a low-thrust transfer trajectory (Morante et al. 2021). 
Approaches to solve the fuel-optimal problem can be generally 
categorized into direct methods, indirect methods and hybrid 
methods. In the direct method, the fuel-optimal problem is 
converted into a parameter optimization problem and then 
nonlinear programming (NLP) can be applied to obtain the 
solution (Hargraves & Paris 1987). Direct methods do not 
require an accurate initial guess because of a large domain of 
convergence. Besides, the optimization variables have physi- 
cally intuitive meanings. However, direct methods require a 
large amount of computation and may not converge to the 
optimal solution. As for an indirect method, the optimization 
problem is usually turned into a two-point boundary-value 
problem (TPBVP) or multi-point boundary-value problem 
(MPBVP) (Haberkorn et al. 2004; Jiang et al. 2012; Chen 
et al. 2018). Indirect methods provide assurances that the first- 
order necessary conditions are satisfied (Morante et al. 2021). 
Nevertheless, indirect methods normally require an appropriate 
initial guess and are usually difficult to solve due to the small 
convergence radius. Hybrid methods are combinations, which 
have the advantages of both methods, where the time histories 
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of costate variables are directly parameterized and the optimal 
solution is obtained through the NLP method (Gao & 
Kluever 2004). 

A multitude of research efforts has been dedicated to 
optimizing low-thrust gravity-assist (LTGA) trajectories (Casa- 
lino et al. 1999; McConaghy et al. 2003; Petropoulos & 
Longuski 2004; Vasile & Campagnola 2011; Jiang et al. 2012; 
Yang et al. 2015). A broad search algorithm with a simplified 
shape-based trajectory model was applied to generate initial 
estimates for Gravity-Assist, Low-Thrust, Local Optimization 
Program (GALLOP) (Petropoulos et al. 2000; Petropoulos & 
Longuski 2004). This algorithm was implemented in the 
software Satellite Tour Design Program for Low-Thrust 
Gravity-Assist (STOUR-LTGA), which automatically searches 
for low-thrust, gravity-assist trajectories. Subsequently, 
STOUR was utilized to conduct a broad search, and the 
resulting LTGA trajectory candidates were optimized using a 
direct method (McConaghy et al. 2003). In some works, 
patched conic trajectories via CP have been taken as initial 
guesses for the optimization of LTGA trajectories, the under- 
lying idea of which was that optimal impulsive transfers can be 
regarded as a limiting case of the fuel-optimal low-thrust 
transfers with no limit on the thrust level (Okutsu et al. 2006; 
Vasile & Campagnola 2011). An automated approach was 
demonstrated to find the number and sequence of swingbys 
(Englander & Conway 2017). In this approach, two nested 
optimization algorithms were combined. The outer loop used a 
genetic algorithm to select the flyby number and sequence 
while the inner loop solved the corresponding sequence of 
interplanetary legs using a direct method. Different from the 
shape-based method and direct method, a practical homotopic 
method which can be applied to solve LTGA trajectories was 
proposed (Jiang et al. 2012). However, this method primarily 
addressed single gravity-assist scenarios, and achieving reliable 
convergence for multi-gravity assist cases remained challen- 
ging due to the large number of unknowns. To deal with this 
issue, low-thrust trajectories with triple swingbys were divided 
into three trajectories and solved by a three-step method (Yang 
et al. 2015). An indirect gradient-based method was also 
developed to automatically find swingbys under the multi-body 
dynamics (Olympio 2008). 

This work is carried out in three steps. First of all, the 
determination of number and sequence of swingbys is based on 
the patched conic trajectories via CP. It is reasonable to apply 
the results of number and sequence for optimization of LTGA 
trajectories (Okutsu et al. 2006; Vasile & Campagnola 2011). 
This method greatly reduces computation time. Second, the 
initial guesses are globally searched in the design space in two- 
body dynamics. The first-order necessary conditions for 
gravity-assist constraints are applied, as derived in the practical 
homotopic method (Jiang et al. 2012). In addition, the 
necessary conditions for optimizing initial and final times are 
derived in this work. An improved cooperative evolutionary 
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algorithm (ICEA) is used (Lei et al. 2013), where the objective 
function combines the propellant mass and constraints with 
penalty factors. Because of the convergence issue for multi- 
gravity assists in the indirect method, the purpose of this step is 
not to achieve precise solutions for the MPBVP, but to obtain 
feasible parameters essential for the subsequent step. Finally, 
the transfer trajectories are further optimized to achieve higher 
accuracy. For transfer from Earth to 2020 XLs, gravitations of 
Venus, Earth, Mars and Jupiter are considered in order to 
ensure the accuracy of the transfer trajectory. However, in 
multi-body dynamics, Euler-Lagrange equations cause sensi- 
tivity issues in the indirect method. The costate variables 
change rapidly over many orders of magnitude during the 
swingby (Olympio 2008). To address this issue during 
swingby, instead of solving the whole transfer trajectory 
directly, a patched-arc model is used. The whole trajectory is 
segmented into interplanetary trajectory legs governed by 
gravitation and thrust of the EP system, and coasting 
planetocentric phases where the spacecraft is only influenced 
by gravitation. The interplanetary trajectory legs remain outside 
the planets’ sphere of influence (SOD, avoiding the rapid 
changes in costate variables during swingby. 

The structure of this paper is as follows. In Section 2, the 
three-step approach for designing the low-thrust multi-gravity 
assist trajectories is introduced, including selecting number and 
sequence of swingbys, calculating basic parameters of the 
trajectories, and optimizing the trajectories in multi-body 
dynamics. Then the three-step approach is applied to design 
the rendezvous mission to 2020 XLs. Details of the rendezvous 
mission are provided in Section 3. The article concludes in 
Section 4. 


2. Methodology 
2.1. Determination of Swingby Number and Sequence 


In the rendezvous mission to an asteroid, the spacecraft 
departs from LEO, followed by a series of selected planet 
swingbys, culminating in the rendezvous with the target 
asteroid. Below is the search method described in detail. In 
this work, the delta-v required to launch from LEO is assumed 
to be supported by the upper stage of a rocket. The departure 
delta-v is defined by Xu et al. (2007) 


Avı = yllv(@o) — velto IP? + vŽo — VLEO; (1) 


where v (tọ) and vg (to) are, respectively, the velocity vectors of 
the spacecraft and Earth in the heliocentric ecliptic reference 


9. and 
FLEO 


frame at the departure epoch, and w% = 


VLEO = a are, respectively, the magnitudes of the 


Earth-escape velocity and actual velocity at the low Earth 
parking orbit with rLeo = 6578.137 km. 


Yang, Xu, & Li 


In order to obtain more feasible solutions, deep-space 
maneuvers (DSMs) are introduced in cases where relative 
velocities at swingby are not exactly equal (Lei et al. 2017; Gao 
et al. 2019). For an impulsive maneuver during swingby, the 
periapsis can be a useful starting point because it is found 
easily and often close to the optimum solution (Gobetz 1963), 
which is demonstrated in Figure 1. 

When there is a DSM at the periapsis, the angles 6, 62 and 0 
can be respectively presented by Gao et al. (2019) 


f 1 
6) = arcsin) —— = }, 
1+ lal? 
Hp 
62 = arcsin E. S g (2) 


fp 
1+ [IP 
Hp 
llv ` Poll 


cos 0 = cos (61 + 62) = —————, 
IIvsllllvsll 


where jp is the gravitational constant of the planet, r, is the 
periapsis radius, and v% and vy, are respectively the hyperbolic 
excess velocities before and after the DSM. The hyperbolic 
excess velocities are calculated by v% = v(tGa) — vp(tca), 
where vp is the velocity vector of the planet in the heliocentric 
ecliptic reference frame and tga is the swingby date. When 
there is no DSM at swingby, the angles have the relationship 
that 6, = 6, = 0/2. According to Equation (2), the value of r, 
can be determined with 


Vy Vas : 1 
jaran Me Nall = arcsin) ——————— 
IIvsllllvssll 1+ Ivl- 
p 
. 1 
+ arcsin)] ——————— |. (3) 


z 
1+ (SIP 
; 


Given a range of r,, Equation (3) is solved using the bisection 
method (Gao et al. 2019). The norms of velocity relative to the 
planet after and before the DSM vę& can be obtained by 
applying the method from Gao et al. (2019) 


vé = hiže + 222. (4) 
lp 


Then the DSM impulse Avps» is 


Tp lp 


ipa = [mer paia [mr +2., 6 


Therefore, once r, is solved, Avpsm can be obtained. 
The rendezvous delta-v is obtained by 


Av: = |v — vat) II, (6) 
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Figure 1. Hyperbolic trajectory during swingby with DSM. 


where v(t) and v,(t,) are respectively the velocity vectors of 
the spacecraft and target asteroid in heliocentric ecliptic 
reference frame at rendezvous epoch. The total delta-v budget 
can be expressed as 


Av = Avı + Av + > Avy (7) 
k=1 


where no is the number of swingbys. Usually, during mission 
planning, constraints should be imposed on the DSM impulse 
due to the engine’s limited thrust capabilities, while constraints 
should also be placed on the departure delta-v due to the upper 
stage’s limit. However, on the one hand, the optimization 
process in this procedure focuses on determining the number 
and sequence of swingbys. On the other hand, the required 
DSM impulse and excessive departure delta-v are feasible with 
EP. Therefore, these constraints are ignored in this process. 


With the total delta-v budget, the impulsive transfer trajectory 
can be optimized using global optimization algorithm ICEA. 
Since the swingby altitude is constrained, the objective 
function is 


no 


Jgequence =Av+ a max (Tp, ink — Vp,k> 0) 
k=1 


+ P; max (t; = Pr ax? 0), (8) 


where p, and p, are penalty factors, and ty is the given 
maximal flight time for the rendezvous mission. For this case, 
the decision variables encompass the departure epoch fo, flight 
time of spacecraft traveling from Earth to the first swingby 
planet ATo and flight times from the k-th swingby planet to the 
next planet or asteroid AT; (k = 1, 2,...,29). Then the swingby 
epochs tga, (k = 1, 2,...,no) and rendezvous epoch tp can be 
determined with those decision variables. Subsequently, the 
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Represent j with 


r 0 1 m- 
J= PiX Np +p x Np ++ pax Np 


Optimize the implusive trajectory of sequence Earth- P, —P, —---—P, 
using ICEA with the objective function J, 


asteroid 


Prg- 


sequence 


Output the value of objective function 
along with the sequence 


Search complete 


Figure 2. Search algorithm of swingby sequence. 


position vectors of the swingby planets are calculated with 
Keplerian elements.* Next, the velocity vectors of the space- 
craft at different epochs, including v (to), v (tp) and v(tga), can 
be evaluated by solving the Lambert problem (Gooding 1990). 
Finally, the total delta-v budget is derived using Equations (1), 
(5), (6) and (7). 

The candidate swingby planets available for swingby are 
denoted as { Po, Pi,...,Px,—1}, where Np is the count of candidate 
swingby planets. The maximum number of swingbys is No. The 
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approach to search the number and sequence of swingbys is 
depicted in Figure 2. In this diagram, for a given swingby number 
No, the total number of swingby sequence combinations amounts 
to (Np)”°. An integer j between 1 and (Np)”° can be represented in 
base-Np with (j)dec = (Dyy—1-+P1P0)N> = Pny—1 X NETH.. 
+p, X NA+ Po X Np. Therefore, every j can be used 
to represent a transfer sequence for Earth- P, — 
Py, — Py, asteroid. Then the ICEA algorithm is applied to 
optimize every transfer trajectory corresponding to the sequence 
Earth- Pp, — Pp, — ... Fp, asteroid. 


Research in Astronomy and Astrophysics, 24:015020 (20pp), 2024 January 


2.2. Optimization of LTGA Trajectory in Two-body 
Dynamics 


2.2.1. Time-free Fuel-optimal Problem for LTGA Trajectory 


In the second step of this work, only the gravitation of the 
Sun is considered. Most derivations in this optimal problem are 
based on Jiang et al. (2012). The motion of the spacecraft is 
governed by the thrust force of the EP system and gravitation 
of the Sun. The equations of motion (EOM) can be expressed 
as 


r= y, 
vp Hs | Tinax U 
r? nm (9) 
— n Tmax U 
m = ME G 
T80 


where r and y respectively denote the position and velocity 
vectors of the spacecraft in the heliocentric ecliptic reference 
frame, r is the Euclidean norm of r, m is the mass of the 
spacecraft, u, is the gravitational constant of the Sun, Tmax and 
Ip are respectively the maximal thrust and specific impulse of 
the EP system, go is the gravitational acceleration of Earth at 
sea level, u € [0, 1] is the thrust ratio and the unit vector a 
signifies the direction of thrust. In the optimization of the 
LTGA trajectory, quantities about length are nondimensiona- 
lized by the astronomical unit au = 149,597,870 km. Quantities 
about time are nondimensionalized by ,/au*/j,, where the 
Sun’s gravitational constant is = p,=1.32712438 x 
10'! km? s~?. Quantities about mass are nondimensionalized 
by the spacecraft’s initial mass mo. 

Generally, the fuel-optimal problem means maximizing the 
final mass, or minimizing the fuel consumption, which is 
defined by 


t 
j= =“ udt, (10) 


where fp and ty are respectively the initial and final times, which 
are both free. The boundary conditions are as follows 


r(to) — re(to) 
Oo = 1 v (to) — [ve(to) + VEw] = 9, (11) 
m(to) — 1 


(12) 


= e — r(t) =0 


TT Tote — walt) = 0" 


where rg is the position vector of Earth, Vg is the hyperbolic 
excess velocity at departure from Earth and r, is the position 
vector of the target asteroid in the heliocentric ecliptic reference 
frame. In addition, for the multi-gravity assist mission, there are 
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also constraints at the swingby points (Jiang et al. 2012) 


r(tca,k) — Fpc(tca,k) 
=| GA,k P,k(ÍGA,k = 0, (k= 1, 2,...,n0), (13) 


lvli — lvl 
Ok = 1— Fp,k/T prak < 0, (14) 


where rp, is the position vector of the k-th swingby planet in 
the heliocentric ecliptic reference frame. To address both 
equality and inequality constraints, numerical multipliers v, €, 
x and « are introduced. Moreover, it is well known that in a 
fuel-optimal problem, the optimal control law has a bang-bang 
structure, where the discontinuity causes the optimal solution to 
converge within a limited domain. To overcome this challenge, 
a numerical continuation method is introduced (Bertrand & 
Epenoy 2002). There have been various approaches to achieve 
this (Haberkorn et al. 2004). A homotopic parameter € is 
utilized to perturb the performance index (Jiang et al. 2012). 
For convenience, define 


no 
G: = X (Xr Yr + KE) +v By E Oo. (15) 
k=1 
The performance index can be modified as 


noT] Ptak 
1= aga fS f [u — ew(1 — u)]dt 


Lp8o k=1 tGa,k 
IGA 
+f [u — eu(1 — u)]dt 
to 


+ “ [u — cu(l — war} + G,, (16) 


'GA.no 


where Ao is a positive numerical factor applied to normalize the 
solution. There is a normalization condition of (Jiang et al. 
2012) 


j + AGIP + SS dixalÊ +) =1, (17) 


k=1 


where A= (A; Ay; Am) is the costate vector. When the 
parameter £ = 1, the performance index is identical to that of 
the energy-optimal problem, which is easier to solve because of 
its control continuity. By the use of Pontryagin’s Maximum 
Principle (PMP), this problem can be transformed into an 
MPBVP. The corresponding Hamiltonian is constructed (Jiang 
et al. 2012) 


H= yt AÈ Bop 4 hata) A Tmax U 
r m 


[u — eu(1 — u)]. (18) 
Tp 8o 


According to PMP, the optimal thrust direction and magni- 
tude, which minimize the Hamiltonian, are determined by 
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Jiang et al. (2012) 


a=, (19) 
IIA 
0 if SF > € 
1 if SF < — 
u=4) "= if SF < a (20) 
— — — if |SF| <e 
2 2e 
where SF is the switching function, which takes the form 
Ip ry m 
SF=1 pSollAvil A (21) 
om Xo 


Substituting the optimal control law presented by Equation (19) 
into the Hamiltonian function, the differential equations 
governing the costate variables, often called Euler-Lagrange 
equations, can be derived 


. ry 3A- r 
Ar = „(> E 5 r) 


r r 
X, = — Ar, (22) 
3 Tmax U 
Xm = =-— llà, ll. 
m 


The transversality conditions for swingbys are (Jiang et al. 
2012) 


Arla) = Arta) — X13, (23) 


Altda) = Xadi + KB, (24) 


Devinrk 


Pinge = Avlak) — Xaple + KA = 0, (25) 


Fp nk 


= + 
Pak = Heltaa d) — Heltta) — Xin, ` VP kCtGA, K) 


! kC = 0, (26) 


+ Xalik — ig)» apaltoan) — —— 
Pmin-* 


where ip = V55,4/ IWS all is the unit vector of v5 ,, and ap, is 
the acceleration vector of the k-th swingby planet. 
Equations (23) and (24) are used to update the position and 
velocity costate vectors. Detailed information for A, B and C 
can be found in Jiang et al. (2012), and they have the forms 


A 


— Tp,k 1 s4 mooo 
~ Teall | Zang Me cos Oiz) — i |. (27) 


B 


— Tpk 1 y= ety et 
~ Weal lam any kk cos Oi) — ip |; (28) 


— 1 1 — of 
c= mf 4sin20/2(1 al mer cos i; ) 


1 si 2 it ig 
+—— (i; — cos ĝi | | 2 
& w) |- vall "Wesel 


I vco.x Il 


` ap k(tGA,k). (29) 
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In addition, the complementary slackness conditions due to the 
inequality constraints are (Jiang et al. 2012) 


kk > 0. (30) 


In Jiang et al. (2012), the initial time and final time are fixed. In 
contrast, the approach in this paper allows for flexibility in both 
initial and final times. For the time-free problem, the first 
differential boundary conditions are given by Hull (2013) 


KkOk = 0, 


ðG: ac. T 
= — = —| — 1 
H (to) Do” A (to) Fen , (31) 
T 
OG: OG: 
-a = Xr 2 = 2 
Ha =— Fe MD | 2 a (32) 


where x ê (r; v, m). Therefore, the transversality conditions 
are obtained as 


Pi s = He(to) — Ar(to) : velto) — Av(to) - de (to) = 9, (33) 
Pos = H: (tp) = Ar(tp) $ Va(tp) = Ay (tp) ag (tr) =0, (34) 
Am(tp) = 0. (35) 


where a, is the acceleration vector of the target asteroid in the 
heliocentric ecliptic reference frame. For the energy-optimal 
problem (¢€=1), the differential equations expressed by 
Equations (9) and (22) are numerically integrated using the 
classic eighth-order Runge-Kutta integrator with a seventh 
order for automatic step-size control, i.e., RKF7(8) (Fehlberg 
1968). However, as the homotopic parameter £ approaches 0, 
the right-hand sides of ordinary differential equations (ODEs) 
experience abrupt variations around switching points. To 
ensure the accuracy of integration with RKF7(8), the bisection 
method is used to detect the switching points (Martinon & 
Gergaud 2010; Zhang et al. 2015). 

For this MPBVP at £= 1, there are 10 + 9ng unknowns in 
total, including initial and final times fọ, t, the positive 
numerical factor Ao, seven initial values of costate vector X(t), 
3D gravity-assist impulse vectors Avg, = V$ — Vo. (kK=1, 
2,...,9), Swingby dates tga x 4D numerical Lagrange multi- 
pliers x, and 1D Lagrange multipliers x. Meanwhile, there are 
10 + 9ng equations as well, consisting of the 6D Equation (12), 
the 1D Equation (33), the 1D Equation (34), the 1D 
Equation (35), the 4D Equation (13) wy, the 1D 
Equation (30) with the application of inequality 
Equation (14) og the 3D Equation (25) ¢).3,, the 1D 
Equation (26) 4, and the normalization condition 
Equation (17). 


2.2.2. Global Search with ICEA Algorithm 


The main purpose of the second step is to obtain feasible 
basic parameters, which can be applied in the optimization of 
LTGA trajectories in the next step. These basic parameters 
include initial and final times fọ and tp respectively, swingby 
dates tca,x, launch vgs, hyperbolic excess velocities v5; and 
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Figure 3. Scheme of patched-arc model. 


periapsis radius r,, at swingbys. Because it would be difficult 
to solve the MPBVP directly, a global optimization method is 
employed to obtain these parameters. Accurate solutions are 
not required in this search. To reduce the computation time, the 
control law based on the energy-optimal problem (e= 1) is 
utilized to search the parameters. With normalization condition 
Equation (17), the search space can be transformed from an 
unbounded space into a hypersphere (Jiang et al. 2012). The 
8+ 5nọ variables Ap, A(to), Xg and sg are normalized with 
7 + 5no angle variables 6. In addition to the variables involved 
in the MPBVP, launch vz, is considered a decision variable in 
the global optimization process. Meanwhile, as the problem is 
constrained, penalty factors are introduced to turn it into an 
unconstrained one. The objective function is expressed as 


IY = Mp + PAIL + Pr max (t¢ = Le ax? 0), (36) 


where y £ (8; to; AT; ATi; Vec; Avex), I S [Eps Ys Res Ges Ph Aml. 
In this search, the penalty factors are p= 1000 and p, = 1000. 


2.3. Design of LTGA Trajectory in Multi-body Dynamics 


2.3.1. Time-fixed Fuel-optimal Problem Without Swingby in 
Multi-body Dynamics 


In the third step, the third-body perturbations of the planets 
in the solar system are considered. The problem without 
swingby is stated in multi-body dynamics here. Given that the 
orbital elements provided by the Minor Planet Center (MPC) 


N OET 4 
Sun 
Table 1 
Heliocentric Orbital Elements of 2020 XL; 
Epoch (MJD) 59649.0 
Semimajor axis, a (au) 1.0006928 
Eccentricity, e 0.3871117 
Inclination, i (Deg) 13.84744 
Longitude of ascending node, 2 (Deg) 153.59764 
Argument of perihelion, w (Deg) 87.98465 
Mean anomaly, M (Deg) 4.6601 


are heliocentric, the EOM are described in the heliocentric 
ecliptic reference frame. Most derivations are identical to 
those in a time-free fuel-optimal problem in two-body 
dynamics, with the exception of the EOM and Euler- 
Lagrange equations. The motion of the spacecraft is governed 
by the thrust force of the EP system, as well as the 
gravitations of Sun, Earth and other planets in the solar 
system. The EOM can be expressed as 


r=y, 
10 
: Hs r—fi Fi Tnax U 
b= -Sr-V i + a, 
bo g (r —rlP r m 
. Tmax U 
m = > 
Lp8o 


(37) 
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Table 2 
Keplerian Elements of the Swingby Planets 
Venus Earth/Moon Barycenter Mars 
Epoch (MJD) 51,544.5 51,544.5 51,544.5 
Semimajor axis, dp (au) 0.72332102 1.00000018 1.52371243 
å (au/Century) —0.00000026 —0.00000003 0.00000097 
Eccentricity, eọ (Rad) 0.00676399 0.01673163 0.09336511 
é (Rad/Century) —0.00005 107 -0.00003661 0.00009 149 
Inclination, Jp (Deg) 3.39777545 —0.00054346 1.85181869 
I (Deg/Century) 0.00043494 —0.01337178 —0.00724757 
Mean longitude, Lo (Deg) 181.97970850 100.46691572 —4.56813164 
L (Deg/Century) 58,517.81560260 35,999.37306329 19,140.29934243 


Longitude of perihelion, Wo (Deg) 131.76755713 102.93005885 —23.91744784 
Ù (Deg /Century) 0.05679648 0.31795260 0.45223625 
Longitude of the ascending node, No (Deg) 76.67261496 —5.11260389 49.71320984 
Q (Deg/Century) —0.27274174 —0.24123856 —0.26852431 
Table 3 
Candidate Swingby Sequences and Performances 
Transfer sequence no Av (km s7!) Transfer sequence no Av (km s7!) 
E-2020 XL; 0 10.33 E-V-E-V-2020 XL; 3 6.31 
E-E-V-E-V-2020 XL; 4 7.00 E-E-V-2020 XL; 2 7.95 
E-V-2020 XL; 1 7.97 E-E-E-V-E-2020 XL; 4 8.27 
E-V-M-V-2020 XL; 3 8.35 E-V-V-2020 XL; 2 8.39 
E-V-E-2020 XL; 2 8.58 E-E-V-M-V-2020 XLs 4 8.76 
E-E-V-E-2020 XL; 3 8.88 
Table 4 boundary conditions are 

Parameters of the Impulsive Transfer Trajectory Corresponding to the Earth— 

Venus—Earth-Venus-2020 XL; Swingby Sequence r(t) = rŠ, v(t) = vě, m(t) = mo. (38) 
Launch date (m/d/y) 1/6/2025 
Venus’ first swingby date (m/d/y) 3/26/2025 r(t7) = rz, v(t% = ve (39) 
Venus’ first swingby altitude (km) 5025.7 
Earth’s swingby date (m/d/y) 2/6/2026 For the optimization of low-thrust trajectory between planets, 
Earth’s swingby altitude (km) 43,125.2 there are no additional intermediate constraints. Applying the 
Venus’ second swingby date (m/d/y) 1/26/2027 : 
Venus? second swingby-altitude Qkin) 215.4 homotopic parameter £, the performance index and corresp- 
Rendezvous date (m/d/y) 12/12/2027 onding Hamiltonian are as follows 
Departure delta-v (km s7!) 4.52 š 
DSM at Venus’ first swingby (km s~') 1.18 x 107° * Tmax 1y 
DSM at Earth’s swingby (km s~!) 4.54 x 1073 I= Ao J: [u — eu(1 — u)]dt, (40) 
DSM at Venus’ second swingby (km s~') 2.19 x 107° eeu) 
Rendezvous delta-v (km s7!) 1.69 
Total delta-v (km s~!) 6.26 He =Ar-vt+A, 
Flight time (days) 1069.42 
Propellant mass ratio (m,/mo) 44.51% 


where u; (i= 1, 2,...,8 for eight planets, i=9 for Pluto and 
i=10 for the Moon) is the gravitational constant of the 
corresponding celestial body, and r; is the position vector of the 
corresponding celestial body. In this step, the states of planets 
are also determined using the JPL ephemeris DE 421. The 


10 
Hs rrr Ti Tmax U 
a kaa ; Q 
| r? A F z) F m 
i=1 
| 
T 


Ao [u — cul — u)]. (41) 


Lp8go Ip 80 


The optimal thrust magnitude and direction are the same 
with those in two-body dynamics, which are determined by 
Equations (19), (20) and (21). Combining the optimal control 
law and the Hamiltonian function, the Euler—Lagrange 
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Table 5 Table 6 


Basic Parameters for Optimization of LTGA Trajectory Corresponding to the 
Earth-Venus—Earth-Venus-2020 XL; Swingby Sequence 


Launch date (m/d/y) 1/3/2025 
Launch vga (kms!) [5.60, 2.37, 1.52] 
Norm of launch vga (km s~!) 6.27 
Venus’ first swingby date (m/d/y) 3/21/2025 
Venus’ first swingby altitude (km) 2969.6 


[12.25, 2.04, 0.22] 
[12.14, — 2.32, 1.96] 
2/3/2026 
33 578.8 
[10.35, — 6.65, 2.32] 
[10.87, — 5.49, 3.11] 
1/21/2027 
845.6 
[13.79, — 1.59, 2.43] 
[11.88, — 2.97, 7.27] 
3/5/2028 
1156.84 


v1 (km s7») 

vs. kms’ 

Earth’s swingby date (m/d/y) 
Earth’s swingby altitude (km) 

Vx. (km s') 

veo (km s7’) 

Venus’ second swingby date (m/d/y) 
Venus’ second swingby altitude (km) 
¥55,3 (km s7) 

v5.3 (kms™’) 

Rendezvous date (m/d/y) 

Flight time (days) 


equations in multi-body dynamics are described as 


. A BAe n à 3A,- — r) 
X = n(ž -=> r)+ n| z ‘ “(rr — ri) |, 
"Ne r° 2 Ilr = ra? Ilr = ra 
X, = -Ar 
. Tmax U 
An = -IA l. 
m 


(42) 


10 


Orbital Elements of Parking Orbit at Launch in Earth-centered Inertial Frame, 
Date and Velocity Vector at the Boundary of SOI in Earth-centered Ecliptic 
Reference Frame 


Epoch (MJD) 60,678.85 


Semimajor axis, a (km) 6578.137 
Eccentricity, e 0 


Inclination, i (Deg) 28.5 
Longitude of ascending node, Q (Deg) 327.904 
Argument of perigee, w (Deg) 284.087 


Mean anomaly, M (Deg) 0 
Departure delta-v, Av, (km s7!) 4.88 

Date leaving Earth’s SOI (m/d/y) 1/5/2025 
Velocity vector leaving Earth’s SOI (km sl) [5.66, 2.40, 1.54] 


In the third step, the initial time and final time are fixed. There 
is only one transversality condition on the mass costate shown 
in Equation (35). The normalization condition is modified as 


fas ACHI = 1. 


For this TPBVP at £ = 1, by the use of normalization condition 
in Equation (43), there are seven unknowns, i.e., the seven 
angle variables transformed from Ap and A(to). There are seven 
constraints including the boundary conditions in 6D in 
Equation (39) and the transversality condition in 
Equation (35). When the homotopic parameter €< 1, the 
positive numerical factor Ay is held constant at the value 


(43) 
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Table 7 
Orbital Elements of Hyperbolic Trajectories at Periapsis, Dates and Velocity Vectors at the Boundary of SOI in Planetocentric Ecliptic Reference Frame 


Venus’ First Swingby 


Earth’s Swingby Venus’ Second Swingby 


Epoch (MJD) 60,755.72 
Semimajor axis, a (km) 2072.93 
Eccentricity, e 5.352 
Inclination, i (Deg) 157.652 
Longitude of ascending node, 2 (Deg) 11.892 
Argument of periapsis, w (Deg) 283.564 


True anomaly, f (Deg) 0 


Date entering SOI (m/d/y) 3/21/2025 
Velocity vector entering SOI (km s75 [12.40, 2.03, 0.23] 
Date leaving SOI (m/d/y) 3/22/2025 


Velocity vector leaving SOI (km sot [12.19, —2.32, 1.97] 


61,074.51 61,426.77 
2523.43 1601.83 
16.834 5.306 
34.353 107.114 
311.230 356.494 
292.636 291.437 
0 0 
2/2/2026 1/21/2027 
[10.42, —6.69, 2.34] [13.96, —1.63, 2.50] 
2/4/2026 1/22/2027 


[10.90, —5.51, 3.12] [11.91, —2.98, 7.29] 


Table 8 
Parameters of LTGA Trajectory Corresponding to the Earth-Venus—Earth— 
Venus-2020 XL; Swingby Sequence 


Departure date (m/d/y) 1/3/2025 
Norm of Vz (1) (km s~') 6.34 
Propellant mass of first leg (kg) 44 
Venus’ first swingby date at periapsis (m/d/y) 3/21/2025 
Norm of velocity vector at the boundary of SOI (km s~') 12.56 
Propellant mass of second leg (kg) 47 
Earth’s swingby date at periapsis (m/d/y) 2/3/2026 
Norm of velocity vector at the boundary of SOI (km s~') 12.60 
Propellant mass of third leg (kg) 22.3 
Venus’ second swingby date at periapsis (m/d/y) 1/21/2027 
Norm of velocity vector at the boundary of SOI (km s~') 14.28 
Propellant mass of fourth leg (kg) 40.8 
Rendezvous date (m/d/y) 3/5/2028 
Flight time (days) 1156.84 
Total propellant mass (kg) 72.2 
Total propellant mass ratio (%) 9.03 


obtained in the energy-optimal problem. There are seven 
unknowns at £< 1, i.e., the initial values of costate vector 
A(to). The normalization condition is not applicable at £ < 1. 
Compared with the MPBVP for low-thrust multi-gravity 
assists, a good initial guess of TPBVP is easier to obtain using 
global optimization. Thus, a numerical method for solving 
nonlinear equations is applied for this problem to attain precise 
transfer trajectories. Here, the ICEA algorithm is utilized to 
search for an initial guess of the energy-optimal problem. The 
objective function is expressed as 
f@* = pA, (44) 
where 3 is the seven angle variables transformed from Ào and 
A(to) according to the normalization condition, and T“ consists 
of the boundary condition in Equation (39) and the transvers- 
ality condition in Equation (35). In this paper, a hybrid solver 
in the GNU Scientific Library (GSL) is used to solve the 
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TPBVP (Powell 1970). The steps of solving the fuel-optimal 
problem can be summarized as (Jiang et al. 2012): 


1. Apply the ICEA algorithm to search for approximate 
solutions of the normalized costate [A,(to); A, (to); Am (to) 
and positive numerical factor A to the energy-optimal 
problem, i.e., € = 1. The integrator is RKF7(8). 

2. Take the approximate normalized solution of step 1 as an 
initial guess to solve the accurate solution of the energy- 
optimal problem. The integrator is RKF7(8), and the 
solver is a hybrid solver in GSL. Go back to step 1 if it 
does not converge. 

3. Fix the positive factor Ao, and decrease the parameter € 
gradually until zero. Solve the new problem using the 
previous solution as an initial guess. The integrator is 
RKF7(8) combined with the bisection method to detect 
the switching points. If it does not converge, shorten the 
decreasing step and try it again. 

4. Output the solutions of the fuel-optimal problem. 


2.3.2. Patched-arc Model for Optimization of LTGA Trajectory 


The patched-arc model is often used for interplanetary 
missions (Melbourne & Sauer 1965; Pierson & Kluever 1994). 
Swingby planets are often seen as massless, and the swingbys 
are considered as instantaneous events (Johnson 1969). How- 
ever, when considering multi-body dynamics, this assumption 
is not adequate. To deal with the rapid changes of costate 
variables during a swingby, the hyperbolic trajectory within the 
planet’s SOI remains unpowered, only subject to gravitation. 
The SOI is applied to place the limit between each low-thrust 
leg. Once the low-thrust legs remain outside the SOI, there are 
no rapid changes in costate variables anymore. The radius of 
SOI of a planet P is calculated by Bate et al. (2020) 


Tsor (45) 


a 

D 
2R 
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where ap is the semimajor axis of the orbit of planet P, and mp 
and M, are the masses of planet P and the Sun respectively. For 
the kth low-thrust leg starting from the boundary of the k — 1th 
celestial body’s SOI to the boundary of the kth celestial body’s 
SOI, the boundary conditions are expressed as 


k k k 
rt) = ret) + Rpt), 
k k k 
v(t) = vp c(t) + Vert), 
k k-1 
maL) = maD), (46) 
k k k 
r(t) = rect) + Rect), 
k k k 
V(t) = vet?) + Veet). (47) 


In the equations above, the superscript (k) denotes the kth 
low-thrust leg, 1‘) is the date when the spacecraft leaves the 
k — Ith celestial body’s SOI, t{ is the date when the spacecraft 
enters the kth celestial body’s SOI, and Rp, and Vp, are 
respectively the position and velocity vectors of the spacecraft 
in the planetocentric ecliptic reference frame. The Oth celestial 
body is Earth. Define m(t®) = ]. Because the no+ Ith 
celestial body is the target asteroid 2020 XLs, there are 
ppt = tf and Reno CTY) = Ve no + er) = 0. 

The position and velocity vectors of a spacecraft at the 
boundary of Earth’s SOI Rp o(t®), Vp o(tQ)) and the corresp- 
onding date 1!) are determined by integrating the state of the 
spacecraft from the parking orbit after the launch supported by 
the upper stage of the rocket at t= tọ. The orbital elements of 
the parking orbit are derived using the launch vg». In addition, 
the orbital elements of the parking orbit are expressed in the 
Earth-centered inertial frame. For the parking orbit, eccentricity 
€park aNd Semimajor axis dpa are both known, and the mean 
anomaly is set as M=O, which is the perigee. There is a 
minimal admissible inclination imin for given vz... When the 
required inclination of parking orbit ip is larger than imin, the 
inclination iS ipa = ip. Otherwise the inclination is 
ipark = imin- The expressions of imin, longitude of ascending 
node Q and argument of perigee w are shown in the Appendix. 
With the orbital elements and departure delta-v given, tQ, 
RET) and VE“ (2) in the Earth-centered inertial frame are 
obtained by integrating the EOM of an unpowered spacecraft 
until ||Re(|| > rsorz. Then the position and velocity vectors 
in the Earth-centered inertial frame are transformed to 
Rp 9 (2), Vp oa) in the Earth-centered ecliptic frame. 

For the swingby, the position and velocity vectors of the 
spacecraft when it leaves or enters the kth planet’s SOI, and the 
corresponding date, are determined by integrating the state of 
the spacecraft forward or backward from the periapsis at 
t= tga.x The orbital elements of the spacecraft at periapsis are 
calculated with v% and periapsis radius r,. Detailed expressions 
of the orbital elements can be found in the Appendix. With the 
orbital elements at the periapsis, 1&7”, Rex (&+)) and 


Ve.(tk*?) are obtained by integrating the EOM of an 


12 


Yang, Xu, & Li 


unpowered spacecraft forward until ||Rpx|| > rsorp. The 
quantities 4, Rp,(t) and Vp,(t) are obtained by 
integrating the EOM of an unpowered spacecraft backward 
until || Rp || = rsorp. The patched-arc model is illustrated in 
Figure 3. The general algorithm of the patched-arc model is 


presented in algorithm 1. 


Algorithm 1. Algorithm of patched-arc model 


Input: initial time tọ, final time ts launch vg.., swingby dates tga,x, hyperbolic 
excess velocities vik at swingbys, periapsis radius rp, (k = 1, 2,...,nọ) 
Output: the fuel-optimal solutions of all low-thrust legs, the whole trajectory 

and the optimal thrust 
Begin 
evaluate ipark, Qparks Wpark in an Earth-centered inertial frame and Av, and 
transform orbital elements into state vector; 
calculate 10, REM), VEQ) in Earth-centered inertial frame with a 


forward integration, and transform state vector into Rpo (td) Vp o( ti) in 


planetocentric ecliptic reference frame; 
for k = 1 to no do 

calculate the orbital elements at periapsis using veg and periapsis 
radius rpk; 

obtain 14), Reeth?) Vorai”), 1, Ret) and Vp) with 
forward and backward integration; 


end 
for k= 1 to nọ + 1 do 
set the boundary conditions with Equations (46) and (47); 
solve the fuel-optimal problem with the given boundary conditions; 
end 
end 


In this section, the design process has been introduced step 
by step. The design process is concluded as: 


1. Selection of the number and sequence of swingbys based 
on Lambert problems. Compute the performances of 
impulsive transfer trajectories for all possible swingby 
sequences according to the process shown in Figure 2. 

2. Search for feasible basic parameters of LTGA trajec- 
tories. Instead of directly solving the MPBVP in an 
indirect method, a global optimization algorithm is used 
to obtain solutions. Those solutions are not accurate 
enough for the MPBVP, but may be helpful to design an 
LTGA trajectory in this three-step approach. Then basic 
parameters are calculated with the solutions. 

3. In order to design a trajectory with better performance, 
the search in step 2 is repeated many times. The 
performances of LTGA trajectories corresponding to all 
solutions in step 2 are briefly estimated with the patched- 
arc approach. However, in the estimation, to save 
computation time, only two-body dynamics is considered 
in low-thrust legs, while planetocentric legs are still 
governed by multi-body dynamics. Finally, the basic 
parameters with the best performance are used in the 
design of the LTGA trajectory in multi-body dynamics. 
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Figure 5. The LTGA trajectory in the heliocentric ecliptic reference frame 
corresponding to the Earth-Venus—Earth-Venus-2020 XL; swingby sequence. 
The orbits of Earth, Venus and the target asteroid 2020 XL; are displayed by lines 
with different colors; the red in the transfer trajectory means the legs where the 
thruster works, and the green in the transfer trajectory represents the planetocentric 
coasting legs. 


3. Results 


In this section, the rendezvous mission to 2020 XL; is designed 
through the three-step approach. First of all, the swingby sequence 
to 2020 XL; is determined. Second, the basic parameters for the 
LTGA trajectory are obtained through a global optimization. 
Finally, a more precise LTGA trajectory for a rendezvous mission 
to 2020 XL; is optimized in multi-body dynamics. 


3.1. Swingby Sequence for Rendezvous to 2020 XLs5 


The orbital elements of the ET 2020 XL; are listed in 
Table 1, which is provided by the MPC.* 


3 https: //www.minorplanetcenter.net/data 
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Given that the synodic periods of Earth relative to Venus and 
Mars are respectively 1.6 yr and 2.1 yr (Takao et al. 2021), in 
the global search for swingby sequence, the departure epoch 
spans from 2024 January 1 to 2028 January 1, a range 
approximately twice the width of synodic periods to account 
for more possible sequences. All the flight times AT, AT, vary 
between 10 and 1200 days. In previous work, it was 
emphasized that Venus and Earth are significant swingbys 
planets for rendezvous to asteroid 2010 TK; (Lei et al. 2017). 
Additionally, considering the large eccentricity of 2020 XLs, 
with its distance from the Sun ranging from 0.6133 to 1.3881 au, 
Mars is also a favorable swingby planet. Therefore, Venus, Earth 
and Mars are designated as potential swingby planets. The 
minimum swingby altitudes for Venus, Earth and Mars are 
uniformly set as 200 km. The Keplerian elements of Venus, 
Earth/Moon barycenter and Mars are listed in Table 2. In this 
search, the Keplerian elements for Earth/Moon barycenter serve 
as an approximation of Earth’s position. The maximum 
allowable total flight time for the rendezvous mission is set as 
tr, = 1200 days. In the optimization of impulsive transfer 
trajectories for all the sequences using ICEA, the subpopulation 
allocated for the particle swarm optimization (PSO) algorithm is 
NP, = 400, while the subpopulation related to the differential 
evolution (DE) algorithm is XNP,= 100. The number of 
generations is Ng=3000. The penalty factors are set as 
Dr= 10° and D= 10°. For this mission, the number of candidate 
swingby planets is Np=3, and the corresponding set is 
{Po, Pi, Po} = {Venus, Earth, Mars}. Given the mission’s total 
flight time constraint of 1200 days, the maximum number of 
swingbys is restricted to No = 5. 

The computations are executed on a desktop personal 
computer with a CPU of 3.7 GHz, and Visual Studio 2019 is 
used. The best 10 results are listed in Table 3, where E denotes 
Earth, M means Mars and V signifies Venus. For those results, 
the constraints on swingby altitude and total flight time are 
satisfied. Moreover, for comparison, the direct transfer from 
Earth to asteroid 2020 XL; is optimized in the same way. The 
result is listed in Table 3 as well. Direct transfer to 2020 XL; 
requires a launch delta-v of 5.93 kms_' and rendezvous delta-v 
of 4.40 kms’. Both of them are too large for the mission. 
According to this result, the multi-gravity assist technique 
greatly reduces the total delta-v for this mission. Mars is not an 
adequate swingby planet for the rendezvous mission to 2020 
XLs. Moreover, using Venus as the final swingby is the best 
option to reach asteroid 2020 XLs;. The swingby sequence 
Earth-Venus—Earth-Venus-2020 XL; is the best for this 
rendezvous mission, where the number of swingbys is 
no = 3. Therefore, in the following parts, the LTGA trajectory 
corresponding to the Earth-Venus—Earth—-Venus-2020 XL; 
sequence will be optimized. Before the optimization of the 
LTGA trajectory, the impulsive transfer trajectory corresp- 
onding to the Earth-Venus—Earth-Venus-2020 XL; sequence 
is further optimized with Ng = 5000, and related parameters are 
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Figure 6. The time histories of the magnitude of thrust and spacecraft’s mass corresponding to the Earth-Venus-Earth-Venus-2020 XL; swingby sequence. 
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Figure 7. The time histories of the direction of thrust in the heliocentric ecliptic reference frame corresponding to the Earth-Venus—Earth-Venus-2020 XL; swingby 


sequence. 


14 


Research in Astronomy and Astrophysics, 24:015020 (20pp), 2024 January 


>o O OQ 
N o © = 


Semi-major axis a, AU 
E 
D 


500 1000 
Flight time, days 


© 


— aak 
(æ) (Sz) 


Inclination i, deg 
ol 


500 1000 
Flight time, days 


oO 


Eccentricity e 


Lon. Ascn. Node 2, deg 


Yang, Xu, & Li 


2 
D 


d 
oa 


oO 
nN 


500 1000 
Flight time, days 


500 
Flight time, days 


1000 


Figure 8. The time histories of orbital elements in heliocentric ecliptic reference frame including semimajor axis, eccentricity, inclination and longitude of ascending 
node corresponding to the Earth-Venus—Earth-Venus-2020 XL; swingby sequence. 


Table 9 
Parameters of Direct Low-thrust Transfer Trajectory from Earth to 2020 XLs 


Departure date (m/d/y) 2/15/2025 
Departure delta-v, Av; (km s7’) 5.00 
Date leaving Earth’s SOI (m/d/y) 2/16/2025 
Norm of Ve (tQ) (km s7') 6.57 
Rendezvous date (m/d/y) 10/9/2027 
Flight time (days) 966.27 
Propellant mass (kg) 166 
Propellant mass ratio (%) 20.8 


listed in Table 4. In Table 4, the propellant mass ratio is 
calculated by Mingotti et al. (2012) 


Av + ad (48) 


Tspe8o 
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where Ipe = 300s is the specific impulse of the CP system, 
and m, is the propellant mass. For the impulsive transfer 
trajectory corresponding to the Earth—-Venus—Earth—Venus- 
2020 XL; sequence, the midcourse maneuvers at swingbys are 
relatively small. The propellant mass ratio is 44.51%. In order 
to save more propellant, low-thrust and multi-gravity assist 
techniques are applied in the following part. The optimal 
solutions of impulsive transfer trajectory are used to provide 
the upper and lower bounds of launch date, swingby dates and 


rendezvous date. 


3.2. Global Search for Basic Parameters 


In this mission, the maximal thrust magnitude of EP is 
Tmax = 60 mN, the specific impulse of EP is Zp = 3000s, the 
gravitational acceleration of [Earth at sea level is 
go = 9.80665 ms and the initial mass of the spacecraft is 
mo = 800 kg. The states of planets are determined using the 
JPL ephemeris DE 421 (Folkner et al. 2009). The number of 
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Table 10 
Position and Velocity Errors of the Low-thrust Trajectory for Different force Models 


Two-body Force 


Gravitations of Sun, Venus, Earth, Mars and Jupiter 


Position error of first leg 2.77 x 107° 
Velocity error of first leg 4.26 x 107° 
Position error of second leg 1.75 x 107° 
Velocity error of second leg 1.62 x 107? 
Position error of third leg 2.96 x 107° 
Velocity error of third leg 4.54 x 107° 
Position error of fourth leg 2.63 x 1077 
Velocity error of fourth leg 4.61 x 10° 


Note. Both position and velocity errors are nondimensionalized. 
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Figure 9. The LTGA trajectory in the synodic coordinate system corresp- 
onding to the Earth-Venus—Earth-Venus-2020 XL; swingby sequence. 


generations Ng = 100,000, NP; = 160, NP, = 40. Referring to 
the optimized parameters of impulsive transfer trajectories, the 
launch date ranges from 1 October 2024 to 1 March 2025, 
flight time of spacecraft traveling from Earth to the first 


16 


4.88 x 1077 
9.69 x 107° 
1.83 x 1075 
1.00 x 1075 
1.59 x 1075 
3.23 x 10> 
1.61 x 107 
2.72 x 10> 


swingby planet ATọ ranges from 50 to 100 days and flight 
times from the kth swingby planet to the next planet or asteroid 
AT, all range from 300 to 500 days. Because the maximal 
gravity-assist impulse is the local circular speed at the minimal 
admissible periapsis radius (Sims 1996), which means 
[Ave kllmax = 4/ Hp k/o ak» all magnitudes of the gravity- 
assist impulse range from 0 to 10 kms '. Considering the limit 
of the upper stage, the departure delta-v Av, is assumed not to 
exceed the critical velocity Vo=5 kms! for a payload of 
800 kg (Lei et al. 2017). Thus, the launch ||vg..|| is no more 
than 6.50 kms! according to Equation (1). The maximum 
total flight time is ty == 1200 days. 

This search process is repeated for 40 times with these 
configurations. There are totally three solutions that can 
provide feasible basic parameters for optimization in the third 
step. The basic parameters are then derived with those 
solutions. The corresponding LTGA trajectories of these 
solutions are briefly estimated using the patched-are approach, 
where only two-body dynamics is considered in low-thrust 
legs. In two-body dynamics, the minimum fuel consumption is 
8.95% with a total flight time of 1157 days. The worst solution 
exhibits a fuel consumption of 10.4% with a total flight time of 
1142 days. The convergence curves of these results are 
depicted in Figure 4. Theoretically, the convergence curve 
should converge to a value smaller than 1 where f(y) = m, and 
T =0. However, because of the critically small convergence 
domain of the multi-gravity assist problem, reaching this value 
through optimization is exceedingly challenging. According to 
the convergence curves, objective functions of these solutions 
all reach 20 ~ 30. A large generation number is necessary for 
the objective function to reach a small value. It can have more 
of a chance to provide basic parameters that lead to converged 
transfer trajectories with better performances. Nevertheless, a 
larger number of generations is much more time-consuming. 
The basic parameters of the LTGA trajectory which has the 
minimum fuel consumption are listed in Table 5. 
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Figure 10. The direct low-thrust transfer trajectory in the synodic coordinate system corresponding to the Earth-2020 XL; sequence. 


3.3. LTGA Trajectory from Earth to 2020 XL5 in Multi- 
body Dynamics 


In this mission, the semimajor axis of the parking orbit is 
Apark = 6578.137 km. The eccentricity is epak = 0. The 
preferred inclination of the parking orbit is ig =28°5. The 
gravitations of the Sun, Venus, Earth, Mars and Jupiter are 
taken into account. With the basic parameters given in the last 
section, the LTGA trajectory is optimized in multi-body 
dynamics with algorithm 1. The orbital elements of the parking 
orbit and departure delta-v are listed in Table 6. The orbital 
elements of the hyperbolic trajectories at periapsis are listed in 
Table 7. Parameters of the LTGA trajectory are presented in 
Table 8. The trajectory of the spacecraft in the heliocentric 
ecliptic reference frame is depicted in Figure 5. Because the 
target asteroid 2020 XL; has a large inclination with respect to 
the ecliptic plane, the final Venus swingby must be performed 
around the date when Venus crosses the ascending or 
descending nodes of 2020 XLs. This opportunity comes every 
112.35 days, as the orbital period of Venus is about 224.7 days. 
This analysis is identical to the trajectory depicted in Figure 5. 
The time histories of the magnitude and direction of thrust and 
spacecraft’s mass are shown in Figures 6 and 7. The time 
histories of orbital elements in the heliocentric ecliptic 
reference frame, including semimajor axis, eccentricity, 
inclination and longitude of ascending node, are displayed in 
Figure 8. As affirmed in Table 8, the total flight time is 1156.84 
days, and the propellant mass ratio is 9.03%. The final mass of 
the spacecraft is 728 kg. For comparison, the direct transfer 


17 


trajectory from Earth to 2020 XLs is optimized with the 
patched-are model as well, whose parameters are presented in 
Table 9. As shown in Table 9, the total flight time is 966.27 
days, which is 190.57 days shorter than the transfer with the 
multi-gravity assist technique. However, the propellant mass 
ratio is 20.8%. It means that the multi-gravity assist technique 
reduces the propellant mass ratio by 11.8%, corresponding to a 
saving of 94 kg in propellant mass. Therefore, the multi-gravity 
assist technique is significant for this mission to save the fuel 
consumption. 

In order to show the influence of force models on the 
trajectory design result, the control laws derived from the two- 
body dynamics model and multi-body dynamics, including the 
gravitations of Sun, Venus, Earth, Mars and Jupiter, are 
respectively applied in the full gravitation system (gravitations 
of Sun, eight planets and the Moon). The position and velocity 
errors at the end of every leg are listed in Table 10. The 
position and velocity errors in Table 10 are both nondimensio- 
nalized. The results show that the control law derived from a 
two-body force model can lead to a large position and velocity 
errors (~10~*) at the end of every low-thrust leg, while the 
control law derived under multi-body dynamics, including 
gravitations of Sun, Venus, Earth, Mars and Jupiter, can keep 
the position and velocity errors around 1075. The inclusion of 
Mercury, planets outside Jupiter and the Moon will greatly 
increase the computation time for optimization of the transfer 
trajectory. Therefore, for the transfer trajectory to 2020 XLs, a 
multi-body dynamics model, including gravitations of Sun, 
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Figure 11. The Sun—Earth Jacobi constant of direct transfer trajectory (above) 
and LTGA trajectory (below). Red color means the parts where the thruster 
works, and green color in the transfer trajectory signifies the planetocentric 
coasting parts. 


Venus, Earth, Mars and Jupiter, not only ensures accuracy, but 
also saves computation time. 

For transfer to the ET asteroid, some characteristics in the 
Sun-Earth synodic coordinate system are significant. The 
trajectory of the spacecraft in the synodic coordinate system is 
depicted in Figure 9. The trajectory of direct transfer in the 
synodic coordinate system is shown in Figure 10 as well. In the 
Circular Restricted Three-Body Problem (CR3BP), the trajec- 
tories that lie along an unstable manifold of an L, orbit 
naturally approach the L4 region (Elliott et al. 2020). In 
Figures 9 and 10, the two transfers both approach 2020 XL; 
through trajectories near the Lı point. The Jacobi value 
represents a constant of motion within the rotating frame 
(Sood & Howell 2019). With the definition of Jacobi constant, 
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a low value of the Jacobi constant is analogous to a higher 
energy for the spacecraft as it moves within the Sun—Earth 
system (Elliott et al. 2020). Although the Jacobi constant is 
only conserved in the CR3BP, and not in any higher-fidelity 
models, it is valuable in preliminary trajectory design activities. 
The spacecraft orbiting around Earth has a large value of Jacobi 
constant (~3). However, the target asteroid 2020 XL; is 
orbiting around the L4 point with a low value of Jacobi constant 
(~2.8). From the viewpoint of energy, the transfer from Earth 
to 2020 XLs is challenging. The time histories of Jacobi 
constant of transfer with and without swingbys are both 
illustrated in Figure 11. For direct transfer, the thruster works 
for most of the flight to reach the desired Jacobi value. In 
contrast, employing multiple swingbys greatly reduces the 
Jacobi constant, requiring the thruster to work for a relatively 
shorter duration. It highlights the significance of the multi- 
gravity assist technique in reducing fuel consumption for this 
rendezvous mission. From the viewpoint of the Jacobi constant, 
Venus’ first swingby reduces most of the Jacobi value (by 
~0.1), and Venus’ second swingby reduces the Jacobi value by 
0.04. Judging from the time histories of orbital elements in 
Figure 8, the semimajor axis, eccentricity and longitude of 
ascending node are mainly changed by Venus’ first swingby, 
while the inclination is greatly changed by Venus’ second 
swingby. 


4. Conclusions 


In this paper, a rendezvous mission to Earth’s second ET 
asteroid 2020 XL; with low-thrust multi-gravity assist 
techniques is proposed. A three-step approach is introduced 
to design the LTGA trajectories in multi-body dynamics. The 
key conclusions of our study are as follows: 


1. The direct impulsive transfer from Earth to 2020 XL; 
requires a total delta-v of 10.33 kms", which is 
challenging for the capacity of current technology. 

2. It is found that the best transfer sequence for a 
rendezvous mission to 2020 XL, is Earth—-Venus— 
Earth-Venus-2020 XL;. For the optimization of impul- 
sive transfer trajectory corresponding to this sequence, 
the total delta-v is 6.26 kms ', with a fuel consumption 
of 44.51%. It highlights the potential of the multi-gravity 
assist technique in reducing both the launch energy and 
rendezvous delta-v for such a rendezvous mission. 
Moreover, Venus takes an important role in the 
rendezvous mission to 2020 XLs. 

3. The direct low-thrust transfer trajectory from Earth to 
2020 XL; requires fuel consumption of 20.8%. This 
aspect emphasizes that the EP system is significant for 
saving fuel consumption in an interplanetary mission. 

4. For the LTGA trajectory corresponding to the Earth- 
Venus—Earth—Venus-2020 XL; sequence, the fuel con- 
sumption is 9.03%, demonstrating the effectiveness of 
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combining EP and multi-gravity assist techniques 
compared to using either individually. Therefore, low- 
thrust multi-gravity assist techniques are recommended in 
the rendezvous mission to 2020 XLs. 


Acknowledgments 


This work was supported by Basic Research Project of China 
(grant No: JCKY2020110C096) and the National Key R&D 
Program of China (grant No: 2020YFC2201202). 


Appendix 


In this appendix, some expressions about the conversion 
between launch v (transformed by vz, from Earth-centered 
ecliptic reference frame to Earth-centered inertial frame) or 
swingby v% and the orbital elements of the parking orbit or 
hyperbolic trajectory are listed. 


For the parking orbit, the minimal admissible inclination is 


, Yoo,21 
Imin = arctan SS i 
2 
Vox T Voy 


where the subscripts x, y, z represent the three components of 
the vector in the reference frame. The longitude of ascending 
node Q is given by 


Woark Xx 
Qpar = arctan] —P***_]. 
— Woark,y 


The quadrant of Qpark is chosen according to the sign of Wpark,x 
and —Wparky. Moreover, W is the unit vector of angular 
momentum, which is defined by 


_ rx 
Ilr x vll 


The x and y components of Wpark are calculated by 


Woark,z = COS Ipark, 
aad 2 

q = Vox + Voo,y? 

qz = 2¥o0,y¥oo,z Wpark,z» 
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The two solutions of Wparx are both applicable. In this work, 
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the solution Wpark,y = 24 is taken. The argument of 
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perigee w is obtained by 


Veo ,x Woark,y = V>0,y Wpark,x 
Jo = arccos : 


[vol [sin Ipark 


f } if vy. >0 
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For the hyperbolic trajectory at swingby, the semimajor axis is 


— _tp 
IIvscIP 

Note that because the solution in Section 3 may not completely 
satisfy the constraints, there is a small error between ||v,,|| and 
I[vtI|, which is like the swingby with DSM illustrated in 
Figure 1. To obtain the hyperbolic trajectory without DSM, the 
orbital elements are evaluated mainly based on the last part of 
the hyperbolic trajectory where the spacecraft is going away 
from the planet, and ||v5|| is used. The eccentricity, inclination 
and longitude of ascending node of the hyperbolic trajectory 
are expressed as 


AGA 


7 
PpP 
eca = 1+ , 
aGA 
iga = arccosWea.-, 
Woa, 
Qoa = arctan| ——— |, 
—Woay 


where the unit vector of angular momentum Wea is calculated 
by 


The quadrant of Qga is chosen according to the sign of WoA, 
and —Wga,y. The argument of periapsis is 


+ + 
— Vox WGA,y + mt) 


Ilvšllsin ica 


f= | 


woa = fi - arceos(——-), if Lane >0 


eGA 


waa = wm -fi — arceos(——-), if vi. <0 
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The true anomaly is zero at the periapsis. 
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